################ Table 3 

sample_tab3=dfbartik[dfbartik$testmonth=="May" & dfbartik$subject=="rla",]

# calculate average district race characteristics
dfrace=unique(sample_tab3[,c("stateabb","leaid","grade","pct_black_2005")]) %>% group_by(leaid) %>% summarize(avgpctblack=mean(pct_black_2005,na.rm=TRUE)) %>% ungroup() 

# get district rankings 
dfrace=dfrace %>% mutate(q_black=ntile(avgpctblack,4)) %>% ungroup()


# merge back in:
sample_tab3_race=merge(sample_tab3, dfrace[,c("leaid","avgpctblack","q_black")],by="leaid")

sample_tab3_race$highblack=ifelse(sample_tab3_race$q_black>2,1,0)


bartik_weathavgsort=as.formula(paste("cs_mn_all~ cs_mn_all_lag_cohort + ", districtcovariates," + ", tempvarsavg ," +",acs," |factor(year) +factor(grade)+factor(leaid) + factor(subject) |(avgpm25_9 | I(avgpm25_9*highblack) ~ ",instname," + I(",instname,"*highblack))|leaid",sep="" ))

resulthighblk=felm(bartik_weathavgsort,sample_tab3_race,weights=sample_tab3_race[[weightval]],na.action=na.omit)

resulthighblk$stage1$iv1fstat



########### INCOME ################# 
dfinc=unique(sample_tab3[,c("stateabb","leaid","grade","pct_FRPM_2005")]) %>% group_by(leaid) %>% summarize(avgpctlowinc=mean(pct_FRPM_2005,na.rm=TRUE)) %>% ungroup() 
dfinc=dfinc %>%  mutate(q_lowinc=ntile(avgpctlowinc,4)) %>% ungroup()


# merge back in:
sample_tab3_inc=merge(sample_tab3, dfinc[,c("leaid","avgpctlowinc","q_lowinc")],by="leaid")

sample_tab3_inc$highFR=ifelse(sample_tab3_inc$q_lowinc>2,1,0)



bartik_weathavgsort=as.formula(paste("cs_mn_all~  cs_mn_all_lag_cohort +", districtcovariates," + ", tempvarsavg," +",acs,"|factor(year) +factor(grade)+factor(leaid) + factor(subject)|(avgpm25_9 | I(avgpm25_9*highFR) ~ ",instname," + I(",instname,"*highFR))|leaid",sep="" ))
resultfr=felm(bartik_weathavgsort,sample_tab3_inc,weights=sample_tab3_inc[[weightval]],na.action=na.omit)

resultfr$stage1$iv1fstat


dfinc_race=merge(unique(sample_tab3_inc[,c("leaid","highFR")]),unique(sample_tab3_race[,c("leaid","highblack")]),by="leaid")

############ special needs 
dfspec=unique(sample_tab3[,c("stateabb","leaid","grade","pct_speced_2005")]) %>% group_by(leaid) %>% summarize(avgpctspec=mean(pct_speced_2005,na.rm=TRUE)) %>% ungroup() 
dfspec=dfspec %>%  mutate(q_spec=ntile(avgpctspec,4)) %>% ungroup()


# merge back in:
sample_tab3_spec=merge(sample_tab3, dfspec[,c("leaid","avgpctspec","q_spec")],by="leaid")

sample_tab3_spec$highspec=ifelse(sample_tab3_spec$q_spec>2,1,0)

dfinc_spec_race=merge(dfinc_race,unique(sample_tab3_spec[,c("leaid","highspec")]),by="leaid")

bartik_weathavgsort=as.formula(paste("cs_mn_all~ cs_mn_all_lag_cohort + ", districtcovariates," + ", tempvarsavg ," +",acs," |factor(year) +factor(grade)+factor(leaid) + factor(subject)|(avgpm25_9 | I(avgpm25_9*highspec) ~",instname," + I(",instname,"*highspec))|leaid",sep="" ))
resultspec=felm(bartik_weathavgsort,sample_tab3_spec,weights=sample_tab3_spec[[weightval]],na.action=na.omit)

resultspec$stage1$iv1fstat

############## esl
dfesl=unique(sample_tab3[,c("stateabb","leaid","grade","pct_ell_2005")]) %>% group_by(leaid) %>% summarize(avgpctesl=mean(pct_ell_2005,na.rm=TRUE)) %>% ungroup() 
dfesl=dfesl %>%  mutate(q_esl=ntile(avgpctesl,4)) %>% ungroup()


# merge back in:
sample_tab3_esl=merge(sample_tab3, dfesl[,c("leaid","avgpctesl","q_esl")],by="leaid")

sample_tab3_esl$highesl=ifelse(sample_tab3_esl$q_esl>2,1,0)

dfinc_spec_race_esl=merge(dfinc_spec_race,unique(sample_tab3_esl[,c("leaid","highesl")]),by="leaid")


bartik_weathavgsort=as.formula(paste("cs_mn_all~ cs_mn_all_lag_cohort + ", districtcovariates, " + ", tempvarsavg ," +",acs,"  |factor(year) +factor(grade)+factor(leaid) + factor(subject) |(avgpm25_9 | I(avgpm25_9*highesl) ~ ",instname," + I(",instname,"*highesl))|leaid",sep="" ))
resultesl=felm(bartik_weathavgsort,sample_tab3_esl,weights=sample_tab3_esl[[weightval]],na.action=na.omit)

resultesl$stage1$iv1fstat
stargazer( resulthighblk,resultfr,resultspec,resultesl, keep="avgpm25_9",digits=4, out="output/Table3.tex")
